#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>
#include<string>

double source(double t){
	return sin(6*M_PI*t);
}

void wave2d(	double xa, double ya, int n, double h, 
		double (*Bxa)( double y, double t ), double (*Bxb)( double y, double t ), 
		double (*Bya)( double x, double t ), double (*Byb)( double x, double t ), 
		int nt, double dt, double (* u_initial)(double x, double y), 
		string file ){
	
	int i,j,k;
	double alpha = dt*dt/h/h;
	int ct = n/2;

	ofstream gnu("wave2d_interference.gnu");
	gnu<<"set term gif animate delay 5\n";
	gnu<<"set output \"wave2d_interference.gif\"\n";
	gnu<<"set xrange [0:10]\n";
	gnu<<"set yrange [0:10]\n";
	gnu<<"set zrange [-0.5:1]\n";
	gnu<<"unset border\n";
	gnu<<"unset xtics\n";
	gnu<<"unset ytics\n";
	gnu<<"unset ztics\n";
	gnu<<"set view 45, 30, 1, 1\n";

	double **u_old = new double *[n];
	for(i=0;i<n;i++) u_old[i] = new double [n];
	double **u = new double *[n];
	for(i=0;i<n;i++) u[i] = new double [n];
	double **u_new = new double *[n];
	for(i=0;i<n;i++) u_new[i] = new double [n];

	// calculate u^0_{i,j}
	for(i=0;i<n;i++){
		u_old[0][i] = Bxa( ya + i*h , 0 );
		u_old[n-1][i] = Bxb( ya + i*h, 0 );
		u_old[i][0] = Bya( xa + i*h, 0 );
		u_old[i][n-1] = Byb( xa + i*h, 0 );
	}
	u_old[ct-15][ct] = source(0);
	u_old[ct+15][ct] = source(0);
	for(i=1;i<n-1;i++){
		for(j=1;j<n-1;j++){
			if((i!=ct-15 || j!=ct) && (i!=ct+15 || j!=ct) )
			u_old[i][j] = u_initial( xa + i*h, ya + i*h );
		}
	}
	// output u^0_{i,j}
	string file_frame = file + "000.txt";
	ofstream fp(file_frame);
	if(!fp){
		cout<<"error: failed to open "<<file_frame<<endl;
		exit(1);
	}
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			fp<< xa + i*h <<"\t"<< ya + j*h<<"\t"<< u_old[i][j]<<endl;
		}
		fp<<endl;
	}
	fp.close();

	gnu<<"splot \""<<file_frame<<"\" with pm3d t sprintf(\"t=%f s\", 0)\n";


	// calculate u^1_{i,j}
	for(i=0;i<n;i++){
		u[0][i] = Bxa( ya + i*h , dt );
		u[n-1][i] = Bxb( ya + i*h, dt );
		u[i][0] = Bya( xa + i*h, dt );
		u[i][n-1] = Byb( xa + i*h, dt );
	}
	u[ct-15][ct] = source(dt);
	u[ct+15][ct] = source(dt);
	for(i=1;i<n-1;i++){
		for(j=1;j<n-1;j++){
			if((i!=ct-15 || j!=ct) && (i!=ct+15 || j!=ct) )
			u[i][j] = (1-2*alpha)*u_old[i][j] + alpha/2*( u_old[i-1][j] + u_old[i+1][j] + u_old[i][j+1] + u_old[i][j-1] );
		}
	}
	
	// output u^1_{i,j}
	file_frame = file + "001.txt";
	fp.open(file_frame);
	if(!fp){
		cout<<"error: failed to open "<<file_frame<<endl;
		exit(1);
	}
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			fp<< xa + i*h <<"\t"<< ya + j*h<<"\t"<< u[i][j]<<endl;
		}
		fp<<endl;
	}
	fp.close();

	gnu<<"splot \""<<file_frame<<"\" with pm3d t sprintf(\"t=%f s\", "<<dt<<")\n";

	for(k=2;k<=nt;k++){
		// calculate u^{k}_{i,j}
		for(i=0;i<n;i++){
			u_new[0][i] = Bxa( ya + i*h , k*dt );
			u_new[n-1][i] = Bxb( ya + i*h, k*dt );
			u_new[i][0] = Bya( xa + i*h, k*dt );
			u_new[i][n-1] = Byb( xa + i*h, k*dt );
		}
		u_new[ct-15][ct] = source(k*dt);
		u_new[ct+15][ct] = source(k*dt);
		for(i=1;i<n-1;i++){
			for(j=1;j<n-1;j++){
				if((i!=ct-15 || j!=ct) && (i!=ct+15 || j!=ct) )
				u_new[i][j] = (2-4*alpha)*u[i][j] - u_old[i][j] + alpha*( u[i-1][j] + u[i+1][j] + u[i][j+1] + u[i][j-1] );
			}
		}
		for(i=0;i<n;i++){
			for(j=0;j<n;j++){
				u_old[i][j] = u[i][j];
				u[i][j] = u_new[i][j];
			}
		}
		// output u^k_{i,j}
		file_frame = file + (char)(48+(k/100)%10) + (char)(48+ (k/10)%10) + (char)(48+ (k%10)) + ".txt";
		fp.open(file_frame);
		if(!fp){
			cout<<"error: failed to open "<<file_frame<<endl;
			exit(1);
		}
		for(i=0;i<n;i++){
			for(j=0;j<n;j++){
				fp<< xa + i*h <<"\t"<< ya + j*h<<"\t"<< u_new[i][j]<<endl;
			}
			fp<<endl;
		}
		fp.close();
		gnu<<"splot \""<<file_frame<<"\" with pm3d t sprintf(\"t=%f s\", "<< k*dt <<")\n";
	}
	gnu<<"set output"<<endl;
	gnu.close();
}

double u_initial(double x, double y){
	return 0;
}

double Bxa(double y,double t){
	return 0;
}
double Bxb(double y, double t){
	return 0;
}
double Bya(double x, double t){
	return 0;
}
double Byb(double x, double t){
	return 0;
}

int main(){

	int i,j,n=300;
	double xa=0, xb=10;
	double ya=xa, yb=xb;
	double h = (xb-xa)/n;
	int nt=300;
	double ta=0, tb=5;
	double dt = (tb-ta)/nt;
	cout<<"dt="<<dt<<endl;

	string file = "wave2d/interference";
	wave2d(	xa, ya, n, h, 
		Bxa, Bxb, Bya, Byb,
		nt, dt, u_initial, file );

	return 0;
}
